Analytical Results for Nontrivial Polydispersity Exponents in Aggregation Models 
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We study a Smoluchowski equation describing a simple mean-field model of particles moving in 
d dimensions and aggregating with conservation of 'mass' s = R D (R is the particle radius). In the 
scaling regime the scaled mass distribution P(s) ~ s~ r , and r can be computed by perturbative 
and non perturbative expansions. A possible application to two-dimensional decaying turbulence is 
briefly discussed. 
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Aggregation phenomena arc widespread in nature. 
They have such an impact on material sciences, chem- 
istry, astrophysics, that a large amount of literature was 
devoted to them |]J]. In such dynamical processes, par- 
ticles or objects as different in geometry and size as 
colloidal particles, galaxies, small molecules, vortices, 
droplets, polymers, can merge to form a new entity 
when they come into close contact or interpenetrate, 
through diffusion (Brownian coagulation |^||), ballistic 
motion (ballistic agglomeration exogenous growth 

(droplets growth and coalescence M) or droplet deposi- 
tion @. 

We are usually interested in the evolution of the statis- 
tical distribution of the 'mass' s, a quantity characteristic 
of each particle, that is conserved in the coalescence pro- 
cess: it can be either the real mass, the volume, the area, 
the electric charge, or any other parameter depending on 
the underlying physics. 

Great advance was achieved when it was proposed || 
and observed both in real experiments and in numerical 
simulations that the distribution N(s,t) exhibits scale 
invariance at large time: 



N(8,t)~S(t)-Pf 



S(t) 



(1) 



where the characteristic mass S(t) diverges as t z when 
t — ► oo, ensuring the oblivion of initial conditions and 
physical cut-off or discreteness, as does the diverging cor- 
relation length of critical phenomena: universality arises 
in dynamics as well, with new universality classes. 

The exponents z and (3 are easily deduced from conser- 
vation laws and physical arguments, but in many cases a 
polydispersity exponent r defined by f(x) ~ x~ T when 
x — > is observed, whose value is nontrivial though uni- 
versal. The prediction of r is still a challenge. 

The earliest tool for tackling the problem is still one of 
the most popular, that is the Smoluchowski equation 0. 
It is a master equation pi for the distribution N(s,t): 



dN{ d S t ,t] = \ J N(s u t)N(s - Sl ,t)K( Sl ,s - Sl ) d Sl 
- N(s, t) J N(si,t)K(s, si) dsi 



where the aggregation kernel K(x,y) is symmetric and 
is characteristic of the physics of the aggregation pro- 
cess on a more or less coarse-grained level. Such kinetic 
equations are usually derived within a mean-field approx- 
imation, but in certain cases it is possible to go beyond 
mean-field limitations investigated by van Dongen |Tl| l 
by including approximately spatial correlation effects in 
the kernel K 

One would like to be able to extract the nontrivial ex- 
ponent t for the specific system from the proper kinetic 
equation. This is not an easy task however: the problem 
was clarified by van Dongen and Ernst |Hj who classified 
the kernels according to their homogeneity and asymp- 
totic behavior: 



K(bx,by) = b x K{x,y) 
K{x,y)~x> 1 y v [y»x) 



(3) 
(4) 



(2) 



For a given physical system, the homogeneity A is easily 
determined using scaling arguments. We consider only 
nongelling systems with A < 1 [Q. For > 0, the expo- 
nent r is trivial and found to be r = 1 + A, whereas for 
/i = 0, r depends on the whole solution / of the scaling 
equation derived from Eq. (S) (see Eq. (||) below) . /i < 
does not lead to any power law behavior but rather to a 
bell-shaped scaling function / [0. In the following, we 
shall focus on the /i = case for which the exponent t has 
been so far only determined numerically by direct simu- 
lation of Smoluchowski equation (not an easy task) [fl5[ , 
by time series fl6| , and of course by direct simulation of 
the physical system described by the considered Smolu- 
chowski equation [|. Considering the ubiquity 
and the importance of the \i = case leading to nontriv- 
ial polydispersity exponents, analytical results would be 
certainly welcome. 

The purpose of this Letter is to show, working with 
a physically relevant simple kernel, that some informa- 
tion about r can be extracted from the kernel itself using 
exact bounds, estimates and expansions around exactly 
solvable kernels. We compare our results to numerical 
studies in the literature [Jl5],[l6| and briefly discuss a possi- 
ble original application to two-dimensional decaying tur- 
bulence. 
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Consider hyperspherical particles in a d-dimensional 
box, of polydisperse radii R with distribution F(R,t), 
evolving the following way: at time t we choose the posi- 
tions of their centers with uniform probability in d-space. 
Then each pair of overlapping spheres of radii R\ and R2 
merges to form a new sphere of radius, 



R = (ijf + R2)* 



(5) 



where D is a parameter with D > d. D can be the actual 
dimension of the spheres, as for instance in the case of 
D = 3 spheres deposited on a d = 2 plane |8| . Once each 
coalescence has been resolved, we have reached t + 1. 

The conserved variable is s — R D and the correspond- 
ing kernel for the equation (|^) is K(x,y) — [x~b + y^) d . 
This kernel has been introduced in many contexts from 
molecular coagulation jl^] to cosmology fllS 17 for spe- 



cific values of d and D, and is one of the most studied 
in the literature p^ , p^ -|20[| although very few analytical 
results are known. This kernel has A = and /1 = 0. 
From the conservation law we get [[l4j (3 = 2 and if we 
plug the scaling form into (|2|) we get z = D/(D — d) for 
d < D and: 

sf'(s) + 2f(s) = f(s) J + °° f( Sl )(s? + s±) d d Sl 

- I fa fifii)f{s - Sl )( S f +(s- Sl )i) d d Sl (6) 

If r > 1 each term of the RHS of Eq. (|h is separately 
divergent and they should be properly grouped [jl3]|lj|]. 

In d = or D — 00, Eq. (0) reduces to the constant 
kernel equation with exact solution fo(x) — 2e~ s and 
foo{s) — 2 1 ~ d e~ s . In the case d = 1, D = 1 an exact 
analytic solution is also known for the time dependent 
equation, with r = 3/2 Eq]. 

Now, for given d and I?, and plugging the expected 
small s behavior /(s) ~ s~ T into Eq. (||), one first gets 
that r< 1 + A = 1 + d/D. Then, matching the behavior 
of both sides of Eq. (0) E3JT1, one finds, 



r = 2 



f{x)x D dx 



(7) 



If a > t — 1 we obtain by multiplying Eq. (g) by x a and 
integrating HMO: 



2(1 -a) / x a f{x) dx = 
J f{x)f{y){x^ + y^) d [x a + y a - (x + y) a ] dxdy (8) 

By studying the large s behavior of Eq. (|^), one 
can show that f(s) decays as c 00 s _ oe _;i with c^ 1 = 
f^ 2 (x~k + (1 — x)~k) d x~T> (1 — x)^idx. 

Exact bounds and estimates - We first show that r > 1 
if d > 1. Suppose r < 1 and consider Eq. (®) with 
a = 0. For d > 1, we have {x~s + y°) d > + 
which leads to J f(x)dx < J f(x)dx J f(x)x^dx, hence 
1 > 2 — r, which is contradictory. Notice that Eq. (0) 



with a = 2 for d — 1 and D — 1 leads to J x 2 f(x)dx = 
2(f x 2 f(x)dx)(J xf(x)dx), and we recover the exact re- 
sult r = 2 — J a;/(x) dx = 3/2 [psj in a very simple way. 

We now introduce an extremely simple method of get- 
ting lower and upper bounds for r. We rely on Eq. (|^) 
valid for a > t — 1. Combining Eq. (0) and (|§|), we get: 



(1-a) 



J g{x,y)dxdy 
J 9{x,y)A(x/y) dxdy 



(9) 



where = (1 + u a - (1 + u)")(l + u^) d /{u a + 

uo) satisfies A(u) = A(l/u) and g(x,y) = {x a y~o + 
x ° y a ) f ( x ) f (y) ■ The ratio in Eq. (||) can then be in- 
terpreted as the inverse of a kind of average of A{x/y) 
with the weight g(x,y). For a given a < g?/-D, we deter- 
mine numerically the lower and upper bounds m a and 
M a of the function A(u). Using Eq. (|TJ), this gives 
2 - (1 - a)/m a < t < 2 - (1 - a)/M a . We then 
choose the best values of a < d/D compatible with 
a > t — 1 leading to the tightest bounds. This allows 
us to greatly improve the exact bounds given in E 19 
for d = 1 and to obtain new such bounds for d > 1. For 
instance for the physically interesting cases (see below) 
(d = 1,D = 2), (d = 1,D = 4) and (d = 2,D = 4) we 
respectively found 1.084 < t < 1.147, 1 < t < 1.075 
(compared to 1 < r < 1.28 and 1 < r < 1.109 in @) 
and 1.25 < t < 1.5. 

It is also possible to obtain good estimates by evalu- 
ating the 'average' in Eq. (^ using a reasonable trial 
weight function g(x, y) instead of the unknown exact 
one. A parameter free choice is obtained by replac- 
ing in the above expression of g(x,y) the exact f(x) by 
x~ T exp(— x) which has the correct leading asymptotic 
for small x (by definition of r) and the expected expo- 
nential large x decay [0^3) . This form is known to be a 
good approximation of the actual f(x) obtained in sim- 
ulations J]15[, and is even obtained in exactly solvable 
models pW . The simplest method is to determine r self- 
consistently from Eq. (||) , for instance with a = d/D, but 
the result depends on the chosen a and may even violate 
exact bounds. A much better and hardly more intricate 
method is to choose a sample of values of a, and min- 
imize an error function measuring the violation of the 
corresponding Eq. (0) jl3| . This method can be system- 
atically improved by allowing for n free 'fitting' parame- 
ters (including r itself) in the trial weight g(x,y). Using 
the simplest \ 2 form for the error function with a trial 
function f(x) = (x~ T + cx~n )e~ x (to take into account 
the exact decay at large x), we obtain jTj|: r « 1.111, 
r w 1.016 and r « 1.431 in the three cases considered 
above. The case d = 1, D > 1 has been numerically in- 
vestigated by means of time series in p6| . The authors 
of |[|] found t « 1.11 - 1.12 for D = 2 (using data in 
the text of |l^]) and to t w 1.06 for D — 4 (as roughly 
extracted from Fig. 3 in |l6[ ) , in fair agreement with our 
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bounds and estimates. We shall find later that our esti- 
mate for d — 2 and d — 4 is in good agreement with a non 
perturbative calculation for r and with a perturbative es- 
timate. Our variational method requires very few CPU 
time and is straightforwardly implemented compared to 
methods used in |il||l6| ]. 

Perturbative expansions for d < 1 - Now we use the 
exactly solvable limits d = and D = oo as a basis for a 
perturbative expansion. 

First, we consider the limit d — > 0. We expand / in 
series of d: f(x) — fo{x) + dfi(x) + 0(d 2 ). A system- 
atic way of expanding r would be to write down a lin- 
ear (self-consistent) differential equation for f\ to solve 
it and plug the result into (Q). This method is used 
in |l3) to compute the next order 0(d 2 ). However, as 
far as the first order is concerned we can get it with- 
out solving for f\. By developing the integral expres- 
sion of r, Eq. (0), we get: r = 2 — J f(x)x~idx — 
-d/D J f (x) \nxdx - dj fi + 0(d 2 ). Now we develop 
both sides of Eq. (||) with a — to get an equation 

for I h : I h = 5 J M x )fa(y)M x ^ + V^)dxdy - 
(I fo)(J h), nence I A = ~ I e~ x ~ y ^x{x^ + y^)dxdy. 
After a simple calculation we get: 

T = 2dj\n(^l+ (^-— ~Y^j du + 0(d 2 ). (10) 

Let us mention that we can generalize this result to any 
homogeneous kernel of the form: (g(x,y)) d , leading to, 
T^2dJ^\ng(l,^)du + 0(d 2 ). 

For D = 1, we get r = 2d + 0(d 2 ), in good agree- 
ment with direct numerical integration of Smoluchowski's 
equation performed by Krivitsky [ fl5[ , who obtained r w 
0.2 for d = 0.1 and r « 0.38 for d = 0.2. This result 
for D = 1 also coincides up to order 0(d) with the best 
inequalities for r that we obtained above (as already ob- 
served in iQ), but not for other values of D p3[ . 

The order 0{d 2 ) requires the computation of f\ which 
satisfies a solvable linear second order equation. This 
cumbersome calculation will be presented in pJ. How- 
ever, in the special case D = 1 it is possible to ob- 
tain explicitly the 0(d 2 ) term by expanding Eq. (||) for 
q = d/D. We obtain, r = 2d + (^ - 4)d 2 + 0(d 3 ) 
|^3| . For d = 0.2, we get r w 0.372 compared to the 
already mentioned t « 0.38, whereas for d = 0.4, we get 
t « 0.686 compared to t k 0.7 found numerically in [fl5|| . 

Now, we perform an expansion in powers of 1/D for 
d < 1, expanding f(x) = foo(x) + ^fi(x) + 2^/2 (a) + 
0(1/ D 3 ). We use exactly the same method: we develop 
Eq. (|J) with a = in powers of 1/D, and plug the result 
in Eq. (0), yielding a vanishing first order term and a 
nontrivial second order term: 

, H TT 2 2- d d(l -d) 1 N 



Once again we were able to obtain a highly nontrivial 
expansion of r without solving for fi and f% themselves, 
although this can also be achieved this way ]l3[|. Note 
that in the limit of large D and small d, Eq. (|10|) and 
(|ll|) coincide up to order 0(d/D 2 ). 

Perturbative estimate for d > 1 - In the case d > 1, we 
have shown that r > 1 and since r < l + d/Z),we see that 
t — » 1 for D — > 00 and finite d > 1. The previous pertur- 
bation is not valid because /1 is non integrable. Never- 
theless we can try to obtain an estimate of r in the follow- 
ing way: we make the ansatz / ~ + c/ s 1+£ e~ s when 
s ->■ 0. We plug it into Eq. fjj) and Eq. (|) for a = d/D, 
and after some algebra |l3| we see that for consistency e 
must be of order 1/D and that c = (1 - 2 1 - d )(d/D - e), 
and eventually that e = k/ D + 0(1/D 2 ) where ac is the 
solution of the nonlinear equation: 

I ^ = j\i +u ^rd V (12) 

This equation always has a solution consistent with the 
exact bound 1 < t < 1 + d/D. For instance in the case 
d = 2, D = 4 we obtain r w 1.462, which is highly con- 
sistent with the previous estimate r « 1.431. Though it 
is still of order 1/D, the obtained perturbative estimate 
depends on the choice of a. a = d/D seems however to 
be the most natural choice. 

In d = 1, c vanishes and we do not learn much. Notice 
that we have shown that all terms of the d < 1 series for 
r in powers of 1/D vanish for d —> 1, as can be seen in 
Eq. ([ll]) for the two leading ones. Thus, the correction 
to r = 1 for large D may be non perturbative in d = 1, 
which would again rule out the estimate r = 1 + 1/2D of 
p(| , which also violates our rigorous inequalities M,13) . 

If we now take the d — > 00 limit in Eq. ( [l2] ) , we obtain, 
r ~ 1 + A — 2~ d \ (A = d/D), a non perturbative behavior 
in d which is to be related to the results below. 

Large d and D - We now present a non perturbative 
calculation in the limit of large d and D, keeping the ratio 
A = d/D fixed. In this limit, the kernel can be written, 

(x* +y^) d = 2 d (xy)?(l + 0(d/D 2 )) (13) 

and surprisingly transforms into the celebrated 'product' 
kernel [HIIj JT0J2G] . Assuming scaling (a still controver- 
sial subject ||15| , [I3[ ), one can easily show that r = 1 + A = 
1 + d/D @ (see also Eq. (||) and @ and the discussion 
below them, as it corresponds to [i — A/2 > 0). We have 
shown that including higher order corrections in power of 
1 /D does not change the value of r such that the correc- 
tion to t = 1 + A is certainly non perturbative. In fact, 
we can estimate this correction by assuming that for fi- 
nite d and D, f{s) ~ c\/ s 1+x ~ £d for s — » 0. Plugging 
this estimate in Eq. (pj) with the limit kernel of Eq. (|l^) , 
we first get Ed ~ 2~ d c\/(l — A). c\ can be determined 
by matching the coefficients of the leading terms in Eq. 
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(0) using the kernel of Eq. ( J13| ) . After a straightforward 
calculation, one gets c\ in the d — > oo limit: 

r = l + A-2 1 -^ 1 , c A = 2(1-A)/ A - 1 (14) 

h = f [«(i - u)r 1_A/2 [m a + (i - u) a - i] 

We thus find a non perturbative (exponentially small) 
correction to r in the large d and large D limit, consis- 
tent with the result obtained above for d > 1 and large 
D. Note that Eq. ( |l4| ) is also consistent with the exact 
result that r — > 1 as D — > oo for finite d > 1 , a result that 
we obtain by setting A = (as I\ diverges). In the case 
d = 2 and D = 4 of interest below, we already gave the 
estimate r « 1.431, whereas Eq. leads to r » 1.428. 

Physical applications of these results to massive parti- 
cle aggregation systems and the generalization to other 
physically relevant kernels (as the one applying to the 
systems described in |j]J|) will be presented elsewhere 
fj"3"f . In this Letter, we would like to present an orig- 
inal application outside of this field of massive particle 
aggregation, namely the dynamics of vortices in two- 
dimensional decaying turbulence. 

Recently, a statistical numerical model has been in- 
troduced |2l|]22| l which describes the dynamics and the 
merging of vortices with the assumption that the typi- 
cal core vorticity u and the total energy E ~ J v 2 d 2 x ~ 
J2i u 2 R\ are conserved (Ri is the radius of the i-th vor- 
tex) throughout the merging processes. This model re- 
produces the main features observed in direct numerical 
simulations (sec [^1 22 for details). For instance, af- 
ter noting that a distribution of vortex radii satisfying 
P(R) ~ R^ 13 is equivalent to a Gaussian energy spec- 
trum E(k) ~ k 13 ^ 6 p2| ; the simulation of this model was 
able to reproduce the fact that starting from a Batche- 
lor spectrum E(k) ~ fc~ 3 (j3 = 3), the system evolves 
systematically to a steeper spectrum E(k) ~ fc -7 with 
7 = 6 — (3 in the range 7 w 3 ~ 5 Q . 

Now, one expects that the collision kernel between two 
vortices is somewhat intermediate between the ballistic 
hard-disk form a ~ [R\ + R2) pTq ] , and the totally un- 
correlated form a ~ (i?i + -R2) 2 (where the probability of 
colliding is proportional to the probability that two ran- 
domly placed vortices overlap, see also below Eq. (j^)) 
fHH - Thus, one can describe approximately the decay 
of vortices due to mergings by means of Eq. (||) with 
1 < d < 2 and D = 4, as two colliding vortices merge 
into a new one with R = (Rf+R^) 1 ^ 4 in order to conserve 
energy and core vorticity. One thus expects a power law 
radius distribution P(R) ~ R~@, with /3 = D(t - 1) + 1 
and t given by our model. We find values of 7 ranging 
from 7 « 3.3 for d = 2 (taking r w 1.43) to 7 w 4.94 (tak- 
ing t « 1.016) for d — 1, in good qualitative agreement 
with observed exponents. As also found in direct simu- 
lations, the actual exponent (and here the value of the 



effective correct d) could depend on the actual initial con- 
ditions (u>, area occupied by the vortices ~ enstrophy). 
Note that the limit Batchelor case 7 = 3, is obtained 
when taking the naive strict upper bound r = 1 + d/D 
with d = 2 and D = 4. 

In conclusion, we have introduced a systematic scheme 
to obtain exact bounds and good estimates for the poly- 
dispersity exponent in aggregation models. We have also 
implemented perturbative and non perturbative expan- 
sions found to be in good agreement with already known 
numerical results when available. Finally, this kind of 
calculations generalizes to other interesting kernels, with 
possible physical applications in the field of droplet depo- 
sition PJl^l or even two-dimensional decaying turbulence 
as briefly mentioned in the present Letter. 

We are very grateful to F. Leyvraz and S. Redner for 
valuable comments on the manuscript. 
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